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Abstract 

Since the introduction of the Black-Scholes model stochastic processes have played an increasingly 
important role in mathematical finance. In many cases prices, volatility and other quantities 
can be modeled using stochastic ordinary diff'erential equations. Available methods for solving 
such equations have until recently been markedly inferior to analogous methods for deterministic 
ordinary differential equations. Recently, a number of methods which employ variable stepsizes to 
control local error have been developed which appear to offer greatly improved speed and accuracy. 
Here we conduct a comparative study of the performance of these algorithms for problems taken 
from the mathematical finance literature. 
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I. INTRODUCTION 



Stochastic processes play an increasingly important role in mathematical finance as ev- 
idenced by the large and growing literature on stochastic volatility models l|, S S 0, S, 



|. Often these theories are expressed in terms of stochastic ordinary 



differential equations (SODEs). Examples include the Cox-Ingersoll-Ross 



White 



Log Ornstein-Uhlenbeck 



fl, Nelson 0] 



HuU- 



Affinelljil and Log-linear models 
of stochastic volatility. Other schemes like ARCH models js], |l3[ use discrete time differ- 
ence equations which can be viewed as approximations to diffusions ^ , and which are often 
favored for computational and other reasons. SODE based models tend to have closer rela- 
tionships to fundamental theory, but have the drawback that analytic solutions are rarely 
known. In general these equations must be solved using numerical approximation schemes. 

Numerical methods for SODEs have a long historyQ] but until recently these algo- 
rithms have not achieved the speed and accuracy characteristic of analogous methods for 
deterministic ordinary differential equations (ODEs)^^. This is partly due to the lack of 
variable-stepsize algorithms which allow for the control of local error, and partly due to a 
lack of sufficiently high order algorithms. The MAPLE Stochastic Package for example, 
fails to include variable-stepsize routines and most methods are of rather low order. Poten- 
tial solutions to both of these problems have been reported in the last few years. Discussions 
of variable stepsize strategies for SODEsjl^, I3| and some basic observations regarding Tay- 
lor expansions for SODEs^^ have led to the emergence of a number of published [2fl| and 
unpublished j2l| variable-stepsize codes. These algorithms also have a number of promising 
additional features such as linear scaling of computational cost with numbers of Wiener 
processes. 

In this manuscript we perform a variety of tests to see whether the algorithms give the 
expected improved performance and accuracy. We refer to the method developed by Wilkie 
and Qetinba§jl^, 0] as SDE9, and the unpublished commercial method 21 1 as ANISE. We 
will not attempt to discuss how these codes work but merely focus on their performance. 
We do not consider other variable-stepsize codes such as the weak method introduced in 
Ref. {2^ since they have restricted domains of applicabihty. 
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The methods SDE9 and ANISE when apphed to an Ito stochastic differential equation 

m 

dXt = aiXt, t)dt + ^ bi{Xt, t)dWu (1) 
1=1 

for an observable Xt with Wiener processes Wu require knowledge of the partial derivatives 
of the solutions, i.e., 

- aiX t) ^ V ^^'^^'^'^ (3) 

All of the problems we consider are formulated with Ito stochastic differential equations and 
we provide these derivatives for each problem. ANISE and SDE9 are more easily applied to 
Stratonovich stochastic differential equations 

m 

dXt = a{Xt, t)dt + H^t, t) o dWit (4) 
1=1 

for which 

BX 

* h{Xt,t) (5) 



dWit 

= aiXt,t). (6) 



dXt 



dt 



Extensions to jump diffusions [12] are also straightforward, but are not considered here. 

Our study shows that both SDE9 and ANISE yield accurate solutions to a wide variety 
of stochastic volatility models. ANISE tends to be about twice as fast as SDE9. For one 
problem we find that ANISE performs hundreds of times faster than SDE9. Both methods 
provide a means of obtaining high accuracy solutions to SODE problems, and may prove to 
be useful quantitative tools for further research in mathematical finance. 

In section II we explore Monte-Carlo convergence of numerically calculated means and 
variances of price and volatilities for seven stochastic volatility models taken from the finance 
literature. Section III examines the accuracies of the algorithms for individual trajectories. 

II. MONTE-CARLO CONVERGENCE TESTS 

Here our goal is to test the accuracy and compare computational performance for the 
ANISE and SDE9 numerical methods for SODEs discussed in the introduction. To do this 
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we compare exact and numerically calculated average quantities like mean price and mean 
volatility for a selection of models from the finance literature. 

For each model we compare known exact average quantities Xt to the numerical av- 
erages X^-PP"^"^- = ^ SjLi ^t''^ computed from individual stochastic evolutions X^-'^ for 
j = 1,...,N, obtained using the SODE methods. We examine convergence to the ex- 
act solution by varying the number of trajectories A^. For each observable Xt we calculate 
the base ten log of the mean relative error, 





Xt- 






max{ Xt 




j^approx. 1 1 



(7) 



and plot this against time. This specific denominator, max{|Xt|, |X^'*'''"~°^'|}, is chosen since 
some observables Xt pass through zero and relative error can therefore blow up. We also 
examine the relative CPU times for the two methods. All calculations were performed on a 
600 MHz Alpha processor with a requested tolerance of 10^^^. 



A. Nelson Model 

"—I 

In the Nelson model|8| the log-price Pt and volatility Vt obey 

dPt = ^tdWit 

dVt = 9{uj-Vt)dt + V2XeVtdW2t, 



where the normally distributed real stochastic differentials are uncorrelated and have dWu = 
and dW^f = dt. Thus, this model has two Wiener processes and two equations. 

The derivatives required by the SODE methods are given in Table lU We employed a time 
step dt = .1 and integrated to 100. As with ODE methods, the intermediate steps taken in 



Xt 


dXt 

dt 


dXt 
dWu 


dXt 
dW2t 


Pt 





^JVt 





Vt 


0{u; - Vt) - xevt 





V2X9 Vt 



TABLE I: Derivatives for Nelson Model. 

ANISE and SDE9 do not necessarily reflect certain aspects of the true solutions such as the 
positivity of Vj. To avoid floating point problems one thus programs ^J\Vt\ rather than 
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The actual solutions returned for Vt will of course satisfy positivity as we will see in section 
III where we explore the accuracy of individual trajectories. 

We explored convergence for the average quantities Pt, var(Fj), Vt, and var(Vt). The 
known exact solutions for these quantities are given by 

Pi 

var(Pi) 
Vt 

var(14) 



We chose parameters 6 = .035, ui = .636, A = .296, and set initial price and volatility to 
Po = -5 and Vo = .029. 

In Fig. [T] we plot the log base ten mean relative error in mean price Pt against time for 
ANISE (in part (a)) and SDE9 (in part (b)). Results are shown for 10^ (dashed curve), 10^ 
(dot-dashed curve) and 10^ (solid curve) trajectories. Convergence with increasing numbers 
of trajectories is good in both cases. For ten million trajectories the averages have a relative 
accuracy of about one part in a thousand. Errors this small are not visible in plots and 
so we do not show the actual solutions. The requirement of millions to tens of millions 
of trajectories for full convergence is typical of systems of SODEs with white noises, since 
Monte Carlo error bounds scale as the inverse square root of the number of trajectories. 

Figure ^ also shows plots of the log base ten mean relative error in variance of price 
var(Pf) against time for ANISE (in part (c)) and SDE9 (in part (d)). Here the convergence 
is slightly better than that for mean price. 

In Fig. 121 we show plots of the log base ten mean relative errors in the mean and variance 
of the volatility Vt. Here the variance has larger error than the mean. Again ANISE and 
SDE9 show similar rates of convergence. 

The CPU times for various numbers of trajectories are shown in Table |H] ANISE takes 
about 7 seconds to compute 1000 trajectories. The SDE9 calculations take about 60% 
longer. This table also shows that both methods scale well with the number of trajectories. 
In other words there are no rare problematic trajectories. 



Po 

UJ — Vn , n, 

ujt + ^ e~^* - 1 



Voe-'' + uj{l 



-et\ 



1 — A 



-26»(l-A)^^ 



2u;{uj - Vo) 
1-2A 



(8) 
(9) 
(10) 

(11) 
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(a) Error in Pf vs. t for ANISE. 



(b) Error in Pt vs. t for SDE9. 




'NelvarP_e5.s 
'NelvarP_e6.s 
'NelvarP_e7.s 



(c) Error in var(Pf) vs. t for ANISE. 



(d) Error in var(Pt) vs. t for SDE9. 



FIG. 1: Error in mean and variance of Pt for Nelson model. 



B. Hull- White Model 



For the Hull-White modelj^, |5|, |6| the log-price Pt and volatility Vt obey SODEs 

dPt = ^/VtdWu 

dVt = {9 - XVt)dt + ^dW2t. 

Thus, we again have two equations and two Wiener processes. 

The derivatives needed by the numerical methods are given in Table ITTTl A time step of 
dt = .1 was used and the equations were integrated to 100. Once again care must be taken 
to program ^J\Vt\ rather than \A^. 

We explored convergence for the average quantities Pt, vai{Pt), Vt, and var(Vt). The 





(a) Error in Vt vs. t for ANISE. 



(b) Error in Vt vs. t for SDE9. 





(c) Error in var(Vt) vs. t for ANISE. 



(d) Error in var(Vt) vs. t for SDE9. 



FIG. 2: Error in mean and variance of Vt for Nelson model. 



exact solutions for these observables are given by 

var(Fi) 



_^(l-e-'^*) + p(At + e-'^*-l) 



r 



(12) 
(13) 

(14) 
(15) 



for this model. 

We set parameters 6 = .03, A = .035, 7 = .0068, and initial conditions Vq = .029 and 
Po = .5. 

In Fig. El we plot the log base ten error in the mean and variance of Pt for ANISE ((a) and 
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^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.69E+01 


O.llE+02 


1.59 


10^ 


0.68E+02 


O.llE+03 


1.62 


10^ 


0.68E+03 


O.llE+04 


1.63 


10^ 


0.68E+04 


O.llE+05 


1.62 


10^ 


0.68E+05 


O.llE+06 


1.62 



TABLE II: CPU times for Nelson Model in seconds. 



Xt 


dXt 

at 


axt 

dWit 


dXt 
dWit 


Pt 





VVt 





Vt 


e-\Vt 
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TABLE III: Derivatives for Hull- White Model. 

(c), respectively) and for SDE9 ((b) and (d), respectively). Dotted curves show the results 
for 10^ trajectories while dot-dashed and solid curves are for 10^ and 10'' trajectories, re- 
spectively. Good convergence with numbers of trajectories is seen. Once again, convergence 
for the variance is slightly better than that for the mean. 

Figure H] shows the log base ten error in the mean and variance of the volatility Vt for 
ANISE ((a) and (c), respectively) and for SDE9 ((b) and (d), respectively). Again good 
convergence to the exact results is observed. The errors in the variance are larger than those 
in the mean. 

Cpu times are compared in Table IIVI for various numbers of trajectories. ANISE takes 
7.5 s to compute 1000 trajectories. Again we observe that SDE9 takes 50 % longer. Once 
again good scaling is obtained for both methods with the number of trajectories. 



C. Cox-IngersoU-Ross Model 

The SODEs for the Cox-Ingersoll-Ross modelj^ are 

dPt = {adt+^tdWu)Pt (16) 



dVt = ti{9 - Vt)dt + aJVt [pdWu + VI - P^dW^t] (17) 
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■HW_mnP_e5.a' 
■HW_mnP_e6.a' 
■HW_mnP_e7.a' 




'HW_mnP_e5.s 
■HW_mnP_e6.s 
'HW_mnP_e7.s 



(a) Error in Pf vs. t for ANISE. 



(b) Error in Pt vs. t for SDE9. 
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■HW_varP_e5.; 
■HW_varP_e6.. 
■HW_varP_e7., 




'HW_varP_e5.s 
'HW_varP_e6.s 
'HW_varP_e7.s 



(c) Error in var(Pf) vs. t for ANISE. 



(d) Error in var(Pt) vs. t for SDE9. 



FIG. 3: Error in mean and variance of Pt for Hull- White model. 



and so we have two equations with two Wiener processes. In this case the volatihty depends 
on both Wiener processes. 

The derivatives needed by the numerical methods are provided in Table |^ A time step 
of dt = .01 was used and the equations were integrated to 10. 



We look for convergence in four observables; mean log-price InP^, variance in log-price 
var(lnPi), mean volatility Vt, and variance of the volatility var(V^). The exact solutions for 
these quantities are given by 



9. 



\jiPt = \nPo + {a--)t + —{Vo-e){e-^'-l] 



(18) 



9 



(J o 



a 



var(lnF,) = [9 - —{p - —)]t + -{Vo - 9){p - —)te 

K 4:K K ZK 
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(c) Error in var(Vt) vs. t for ANISE. (d) Error in var(Vt) vs. t for SDE9. 

FIG. 4: Error in mean and variance of Vt for Hull- White model. 



(19) 
(20) 

e""*)' (21) 

The parameters were set to a = .1, k = .29368, 9 = .07935, a = .11425, p = —.2 and price 
and volatility was set to initial values Pq = 1 and Vq = .1. 

In Fig. El we show the log base ten error in InPt and var(lnPt) plotted against time for 
ANISE ((a) and (c)) and SDE9 ((b) and (d)). Dashed, dot-dashed and solid curves represent 
errors for runs of 10^, 10^ and 10^ trajectories, respectively. Good convergence is seen for 
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{^[(H-^)(p-^)-% 

2 

Kt 1 \2 



Vt 



= He-"* + ^^(1 - 
Voa' 



If 



-Kt 



e 

-2K.t 



K 



+ ^(1 

2k ^ 



^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.76E+01 


O.llE+02 


1.47 


10^ 


0.76E+02 


O.llE+03 


1.44 


10^ 


0.76E+03 


O.llE+04 


1.46 


10^ 


0.76E+04 


O.llE+05 


1.46 


10^ 


0.76E+05 


O.llE+06 


1.46 



TABLE IV: CPU times for Hull- White in seconds. 





dXt 

at 




dXt dXt 

dWit dW2t 


Pt 






PtVVt 


Vt 


t^{e - Vt) 







TABLE V: Derivatives for Cox-Ingersoll-Ross Model. 



both methods at all times, although errors in the variance are larger than those in the mean. 

Figure El plots errors in Vt and var(Vt). Again we see excellent convergence in both cases. 
Errors in the variance are bigger than those in the mean. 

The cpu times for various numbers of trajectories are shown in Table IVII ANISE takes 
less than 5 s to compute 1000 trajectories. The ratio of cpu time for SDE9 to that of ANISE 
is now a much larger 2.5. This relative slowing down of SDE9 is probably caused by the fact 
that Vt now depends on two Wiener processes. Once again the ratio is independent of the 
number of trajectories indicating that both methods handle all trajectories equally well. 



# Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.48E+01 


0.12E+02 


2.45 


10^ 


0.48E+02 


0.12E+03 


2.46 


10^ 


0.48E+03 


0.12E+04 


2.46 


10*^ 


0.48E+04 


0.12E+05 


2.47 


10^ 


0.48E+05 


0.12E+06 


2.47 



TABLE VI: CPU times for Cox-Ingersoll-Ross model in seconds. 
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'CIRmnlnP_e5.a' 
ClRmnlnP efi.a' 
CIRmninP e7.a 





'GIRmnlnP_e5.s 
GIRmnlnP Rfi.s 
ClRmnlnP e7.s 




(a) Error in In Pt vs. t for ANISE. 



(b) Error in In Pt vs. t for SDE9. 




■CIRvarlnP_e5.s 
■CIRvarlnP_e6.s 
■CIRvarlnP_e7.s 





(c) Error in var(lnPf) vs. t for ANISE. 



(d) Error in var(lnPf) vs. t for SDE9. 



FIG. 5: Error in mean and variance of InPt for Cox-Ingersoll-Ross model. 



D. Log-Ornstein-Uhlenbeck Model 



The fourth example is the Log Ornstein-Uhlenbeck model|7| for price Pt and volatihty 
Vf. In this model 



dPt = iadt + e^'dWu)Pt 

dVt = (a - hVt)dt + ]^[pdWu + v^l - pHW2t] 

and so we again have two equations and two Wiener processes. The volatility depends on 
both Wiener processes. 

The derivatives required by the SODE methods are given in Table IVIII A time step of 
10~^ was used and the equations were integrated to 0.1. 
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'CIRmnV_e5.a' 
'OIRmnV eG.a' 
ClRmnV e7.a' 



I 




(a) Error in Vt vs. t for ANISE. 



(b) Error in Vt vs. t for SDE9. 



1' 







i4 ^ f 



'CIRvarV e5.a' 
'CIRvar\re6.a' 
'CIRvarV_e7.a' 




(c) Error in Yax{Vt) vs. t for ANISE. 



(d) Error in Yav{Vt) vs. i for SDE9. 



FIG. 6: Error in mean and variance of Vt for Cox-Ingersoll-Ross model. 



Xt 


dXt 

at 


dXt 
dWit 


axf, 

dW2t 


Pt 


aPt-\Pte^^{\p + e^^) 


e^^Pt 





Vt 


a + bVt 







TABLE VII: Derivatives for Log Ornstein-Uhlenbeck Model. 



We look for convergence in three quantities; mean log-price InPt, mean volatility Vj, and 
variance in volatility var(Vt). Exact solutions for these observables are given by 

hT^ = \nP^ + at-]^f^dt' exp{2[K)e"''*' + ^(l-e-''*')] + ^(l-e-'''*')} (22) 
Vt = He-'-^ + ^ll-e-^*) (23) 
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(24) 



The solution for In Pf was obtained using a variable-stepsize Runge-Kutta code for ODEsQj. 
Parameters were set to a = 70, 6 = 100, p = .2 and initial conditions Pq = .5 and Vq = .029 
were used. 




■LOUrrinlnP_e5.a 
■LOUrr^nlnP_e6.a 
•LOUrr^nlnP_e7.a 




■LOUmnlnP_e5.s' 
■LOUmnlnP_e6.s' 
■LOUmnlrP_e7.s 



(a) Error in Pt vs. t for ANISE. 



(b) Error in Pt vs. t for SDE9. 



FIG. 7: Error in mean of InP^ for Log-Ornstein-Uhlenbeck model. 



In Fig. 121 we plot the log base ten relative error in In Pt for ANISE in (a) and SDE9 in 
(b). Errors are shown for averages over 10^ (dashed curve), 10® (dot-dashed curve) and 10^ 
(solid curve) trajectories. In all cases a spike in error is seen near the time t = .01 where the 
exact In Pt passes through zero. The absolute error is small and so the spike in relative error 
indicated in the plots is essentially fictitious and convergence is in fact good at all times for 
both methods. 

Figure IHI plots the mean relative error in the mean and variance of the volatility for ANISE 
((a) and (c), respectively) and SDE9 ((b) and (d), respectively). The error in the variance is 
larger than that in the mean. Here we see good convergence for both methods at all times. 

The cpu times for various numbers of trajectories are shown in Table IVnTl ANISE takes 
about 5 s to compute 1000 trajectories. The ratio of cpu time of SDE9 to ANISE is again 
larger than two. The ratio is roughly independent of the number of trajectories. 
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.11 




3UmnV_e5.a' 
5UmnV_e6.a' 
3UmnV_e7.a' 




I Hi 



0.02 0.04 




(a) Error in Vt vs. t for ANISE. 



(b) Error in Vt vs. t for SDE9. 



■LOUvarV_e5.a' 
■LOUvarV_e6.a' 
■L0UvarV_e7.a' 




(c) Error in var(Vt) vs. t for ANISE. 



(d) Error in var(Vt) vs. t for SDE9. 



FIG. 8: Error in mean and variance of Vt for Log-Ornstein-Uhlenbeck model. 



E. Affine Two Volatility Factor Model 

The equations for the log-price and volatilities of an affine two volatility model 

are 
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dPt = fidt + ^^0 + ^iVu + ^2V2t dW^t (25) 
dVu = (ttio + OiiiVu)dt + dWu (26) 
dV2t = (^20 + a2iV2t)dt + dW2t (27) 

and so we now have three equations and three Wiener processes. 

The derivatives required by the numerical methods are given in Table HXl A time step of 
dt = .001 was employed and the equations were integrated to 0.5. 

15 



^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.53E+01 


0.13E+02 


2.46 


10^ 


0.52E+02 


0.12E+03 


2.38 


10^ 


0.52E+03 


0.12E+04 


2.29 


10^ 


0.52E+04 


0.12E+05 


2.28 


10^ 


0.52E+05 


0.12E+06 


2.29 



TABLE VIII: CPU times for Log-Ornstein-Uhlenbeck model in seconds. 



Xt 


dXt 

at 




dXt 
dWit 


dXt 
dW2t 


dXt 
dWit 


Pt 












VCO - ilVlt - i2V2t 


Vit 


aio + aiiVit - 


1 

4 










V2t 


"20 + "21^2* - 


1 
4 











TABLE IX: Derivatives for Affine Two Volatility Factor Model. 
We examined quantities Pt, var(P(), Vu, and var(Vit) which have exact solutions given 



by 



Po + lit 



var(Pj) = [^0 H (^lo H )(e - IJ 



vax{Vu) 



ail oi2i 

+ ^(V^2o + ^)(e"-* 
"21 a2i 

aii 



ail 



ail 



'lot 

Qio 
2ali 



1)2 + Y}^f^2^..t 



(28) 

(29) 
(30) 

(31) 



an 

The parameters were set as = .01, = .1258, ^2 = -0344, /i = .02, aio = .2894, 
ail = 17.4321, 020 = -0602, 021 = 13.6036, Pq = 1, = .2, and V20 = .2. Note that on 
average Vu increases exponentially with a large exponent, and so the noises in the equation 
for the price Pt are strongly weighted. 

In Fig. Elwe show convergence via the log base ten relative error in the mean and variance 
of the price for ANISE ((a) and (c), respectively) and SDE9 (b) and (d), respectively). For 
ANISE the plots show three curves corresponding to runs with averages over 10^ (dashed 
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I ■AI1_mnH_e5.a' 

■Att_mnP_e6.a' 
I ■Aff_mnP_e7.a' 

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



(a) Error in Pf vs. t for ANISE. 



(b) Error in Pt vs. t for SDE9. 



■AII_varP_e5.a' 
■A1l_varP_e6.a 
■Aff_\/arP_e7.a' 




0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 













Aff varP e5.s' 
'Aff_varP_e6.s' 




A 










I 
















i s 




T 

s 




1 ' 













0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



(c) Error in var(Pf) vs. t for ANISE. 



(d) Error in var(Pt) vs. t for SDE9. 



FIG. 9: Error in mean and variance of Pt for Affine model. 



curve), 10^ (dot-dashed curve), and 10'' (solid curve) trajectories. For SDE9 the plots show 
just two curves corresponding to runs with averages over 10^ (dashed curve) and 10^ (dot- 
dashed curve) trajectories. In both cases good convergence is observed toward the exact 
solution. The error in the variance is larger than that in the mean. As we discuss below the 
relative cpu time for SDE9 is much larger than for previous problems. Indeed, the run with 
10^ trajectories did not finish and so does not appear in the figures. 

Figure ITUl plots the errors in Vu and var(Vit) against time for ANISE ((a) and (c), respec- 
tively) and SDE9 (b) and (d), respectively) for the same numbers of trajectories as in the 
previous figure. Error in the variance is larger than that in the mean. Good convergence is 
again observed for both methods. 
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■Aft_mnV_e5.a' ! 

■AII_mnV_e6.a' ■ ' ■Aff_mnV_e5.s' 

^ ^ ^ ^ ^ ^ •Aft_mnV_e7.a'^ | [ ^ j !::■ ■ : : ; ■ ^ 'Aff_mnV_e6.s'^ 

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 " 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



(a) Error in Vu vs. t for ANISE. 



(b) Error in Vu vs. t for SDE9. 




0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



(c) Error in var(Vit) vs. t for ANISE. 



(d) Error in var(T4t) vs. t for SDE9. 



FIG. 10: Error in mean and variance of Vu for Affine model. 



The cpu times for various numbers of trajectories are shown in Table |3 ANISE takes 
about 3.5 s to compute 1000 trajectories. In spite of the fact that good convergence was 
observed for the SDE9 method its computation times show a large jump from the 10'^ 
calculation to the 10^ calculation. For 10^ and 10^ ANISE is several hundred times faster 
than SDE9. Rare trajectories with difficult stochastic paths appear to be responsible for the 
poor performance of SDE9. 
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^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.34E+01 


O.lOE+02 


2.96 


10^ 


0.34E+02 


0.92E+02 


2.74 


10^ 


0.34E+03 


0.88E+05 


261.61 


10^ 


0.34E+04 


0.13E+07 


386.63 


10^ 


0.34E+05 


NA 


NA 



TABLE X: CPU times for Affine model in seconds. 
F. Log Linear Two Volatility Factor Model Without Feedback 

The price and volatilities obevjl^ 



dPt = {{aiQ + ai2V2t)dt + e 
dV2t = a22V2tdt + dW2t 
dV3t = a33V3tdt + dW3t 



l-tlj!3dW,t + A3dW3t]}Pt 



(32) 
(33) 
(34) 



and so we have three equations and three Wiener processes. 

The derivatives needed by the numerical methods are given in Table IXll A time step of 
10~^ was used and the equations were integrated to 0.1. 



Xt 


dXt 

at 


dXt 


dXt 
dW2t 


dXt 
dWst 


Pt 


[aio + ai2^2i - - \f3i3i^izFt]Pt 


VI - (Vis)' FtPt 





i^isFtPt 


V2t 







1 





V3t 










1 



TABLE XI: Derivatives for Log Linear model without feedback. Here Ft = e^^o+^i^Vat 



We calculated InPf, V2t, Vst, and var(V3f), some of which have known exact solutions 



V2t 

var(V3t) 



T/ „a22t 

V2ne 



'30t 
1 



j2033t 



2a 



33 



(35) 
(36) 
(37) 
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Once again we had to solve an ODE 

= a,o + ai2X^2oe"-* - - exp{2p^o + 2/?i3l^3oe"^^* + ^(e^"-* _ i)} (38) 
at I 0^33 

numerically to find InP^. This was again accomplished using a Runge-Kutta algorithm for 
ODEs0- 

The parameters were set as ctio = .0337, ayi = .4820, 0^22 = 1.0043, 0:33 = .0291, 
Ao = 1.0294, /?i3 = .0261, ^13 = .3285, Po = 1, V20 = .1, Vso = .05. 





(c) V2t VS. t for ANISE (d) V2t vs. t for SDE9 

FIG. 11: Means of InPj and V2t for Log Linear model without feedback 



In Fig. ^2 we plot the log base ten relative error in InPt and V2t for ANISE ((a) and (c), 
respectively) and SDE9 ((b) and (d), respectively). In all cases the dashed curve represents 
an average over 10^ trajectories while the dot-dashed and solid curves are for 10^ and 10^ 
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trajectories, respectively. Good convergence is seen in all cases except near t = for InP^. 



The exact solution for InPj vanishes at f = for our initial condition, and poor relative 
accuracy is seen as a consequence. In fact the absolute accuracy is good at all times for 10^ 
trajectories. 



\ 



■LL2Vn_r 
■LL2Vn_mhU3_e 
■LL2Vn_m iU3_e7.a' 



Wy*, 



hi 




(a) Vat vs. t for ANISE 



(b) Vat vs. t for SDE9 





(c) var(y3t) vs. t for ANISE 



(d) var(V3t) vs. t for SDE9 



FIG. 12: Error in mean and variance of V^t for Log Linear model without feedback 



Figure IT^ plots the errors in V^t and var(V3t) against time for ANISE and SDE9. Once 
again, good convergence is observed. Errors in the mean and variance are comparable. 

The cpu times for various numbers of trajectories are shown in Table IXTTl ANISE takes 
6.5 s to compute 1000 trajectories. Once again the ratio of cpu time for SDE9 to that 
of ANISE is a little greater than two and this number is independent of the number of 
trajectories. 
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^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.65E+01 


0.14E+02 


2.11 


10^ 


0.64E+02 


0.14E+03 


2.13 


10^ 


0.64E+03 


0.14E+04 


2.18 


10^ 


0.64E+04 


0.14E+05 


2.18 


10^ 


0.64E+05 


0.14E+06 


2.17 



TABLE XII: CPU times for Log Linear model without feedback in seconds. 
G. Log Linear Two Volatility Factor Model With Feedback 

The equations for this modeljl^ are 

dPt = {{aio + ai2V2t)dt 

+ e'^^°+'^^«^^*+^^*^'*' [^1 - ijf^ - iPl^dWu + 'ipndW^t + ipidWuWPt (39) 

dV2t = a22V2tdt + dW2t (40) 

dVst = aMtdt + (1 + f3ssVst)dWst (41) 

dVu = a^^Vudt + (1 + l3iiVu)dWu. (42) 

In this case we have four equations and four Wiener processes. 

The derivatives required by the numerical methods are given in Table IXIIII A time step 
of dt = 10"'* was employed and the equations were integrated to 0.1. 



Xt 


dXt 

at 


axt 

dWu 


axt 

dWit 


axt 
awst 


axt 


Pt 


{aw + ai2V2t - \Fl 


pFtPt 





i^isFtPt 


i'uFtPt 




-^[V'l3/9l3(l + PMt) + \^IaPia{1 + (3AAVAt)]Ft}Pt 










V2t 


022^2* 





1 










"33^3* - 1/333(1 + /333V3t) 








l+p33V3t 





Vu 


044^4* - ^/344(1 + I^AAVAt) 











1 + pAAVAt 



TABLE XIII: Derivatives for Log Linear model with feedback. Here Ft = e'^io+'^i3V3t+/3i4V4t ^nd 

P= \/'^-^13-^lA- 
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We examined quantities InP^, V2t, V^t, and var(V3() some of which have exact solutions 

(43) 
(44) 

"33 + Ph 

+ , (e(^"33+/3|3)^-l). (45) 





— 1^206 




= ^3oe"^^* 







We obtained In Pt numerically by solving the ordinary differential equation 
(MPt ..... e^/^io- 



dt 



aw + ai2l/2oe"^^* - — -e^Aa^a* e2/5i4^« (46) 



using a variable-stepsize Runge-Kutta scheme 



15l |. The averages e^^"^»' for i = 3,4 were 



obtained from the moments (V^t)" using e^^'* = I]J^o ^(^i*)" (numerically truncated after 
n = 20) and iteration using 



Vu = y.oe"''* (47) 

(T^ = (\/,o)2e(2"-+/5?.)* + i^!il:^(e(2a..+/3?,)t _ g^,.*) _^ ^ (g(2«..+/3a)f _ ^) ^4g) 

Q^M ~l~ Pa ^cijj -|- /jjj 

(T^ = (l/.o)"e("">»+'^^?-)* + n(r2 - l)l3u f dt' (X¥)^ e("°''+'^^^«)(*-*') 



t 



+ ^^^^^ /* dt' (l^;)^ e(""-+^'^a)(*-*'), for n = 3, 4, . . . (49) 
2 JO 

which are also readily obtained using an ODE code. 

The parameters were set as aio = .0279, ai2 = .7281, 022 = 5.9997, 033 = .1227, 
au = 8.2119, (3io = .0486, (3i3 = -0695, = -3130, (333 = -3672, = -3655, V13 = -1077, 
ipu = .0564, with initial conditions Pq = 1, V20 = -1, V30 = -05, V40 = .2. 



In Fig. El we plot the log base ten relative accuracy of In Pt and V2t against time for 
ANISE ((a) and (c), respectively) and SDE9 ((b) and (d), respectively). The dashed curve 
represents an average over 10^ trajectories, while the dot-dashed and solid curves represent 
calculations with 10^ and 10^ trajectories, respectively. In all cases convergence is good 



except for InP^ in the vicinity of zero where the exact solution vanishes and the relative 
accuracy becomes poorly defined. 

Figure El plots the error in V^t and var(V3t) for ANISE and SDE9. Good convergence is 
observed in all cases. Errors in the mean are greater than those in the variance. 

The cpu times for various numbers of trajectories are given in Table IXTVl ANISE takes 
7.5 s to compute 1000 trajectories. Once again ANISE is about twice as fast as SDE9. 
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■LL2VmnlnP_e5.a' 
'LL2VmnlnP_e6.a' - 
'l-l-2VmnlnP_e7.a' - 

0.08 

(a) InPt vs. t for ANISE 

' ' ■LL2VmnU2_e5.a' ■ 

■LL2VmnU2_e6.a' 
'LL2VmnU2_e7.a' ■ 




'LL2VmnlnP_e5.s' 
'LL2VmnlnP_e6.s' ■ 
'LL2VmnlnP_e7.s' 



(b) \nPt vs. t for SDE9 




(c) V2t vs. t for ANISE 



(d) V2t vs. t for SDE9 



FIG. 13: Means of In Pt and V2t for Log Linear model with feedback 



III. ACCURACY FOR INDIVIDUAL TRAJECTORIES 



Here we again request a relative accuracy of 10~^^ and determine what accuracy is in 
fact obtained on average for individual trajectories. While it is unlikely that results of this 
high precision would be required in actual financial applications, it is worth exploring this 
issue for a few problems where exact solutions of the SODEs are known. We find that the 
calculations are not very sensitive to the requested tolerance, and accuracies of 10~^^ are 
sometimes achieved even when the requested tolerance is only 10~^. The calculations are 
also insensitive to the stepsize. 

For each realization of the observable we thus calculate an exact solution Xf and an 
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(a) Vat vs. t for ANISE 



(b) V3t vs. t for SDE9 




'LL2VvarU3_e5.a' ■ 
■LL2VvarU3_e6.a' 
■LL2VvarU3_e7.a' - 



'LL2VvariJ3_e5.s' 
'LL2VvarU3_e6.s' - 
'LL2VvarU3_e7.s' 



0.02 0.04 



(c) var(F3t) vs. t for ANISE 



(d) var(F3t) vs. t for SDE9 



FIG. 14: Error in mean and variance of V^t for Log Linear model with feedback 



approximate solution x^^"^' from which we compute the log base ten relative error 

1^ ^approa;.| 

^""^^^^maxilX^uk"-!}^' 
We plot the average of this quantity against time t for each model. The exact solutions 
involve some difficult integrals which are also computed using the numerical method, so our 
tests are essentially self-consistency checks. 

For both models we have requested large time steps and integrated to very long times 
in order to make the calculation somewhat challenging. The errors shown are computed as 
time averages over short intervals since there are high frequency fluctuations in the data 
which make identification of the line types in the figures difficult. 
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^ Trajectories 


ANISE CPU Time 


SDE9 CPU Time 


CPU Time Ratio SDE9/ ANISE 


10^ 


0.74E+01 


0.15E+02 


2.04 


10^ 


0.75E+02 


0.15E+03 


2.06 


10^ 


0.74E+03 


0.16E+04 


2.08 


10^ 


0.75E+04 


0.15E+05 


2.08 


10^ 


0.75E+05 


0.15E+06 


2.08 



TABLE XIV: CPU times for Log Linear model with feedback in seconds. 



A. Vasicek interest rate model 

The SODE for this model is 

dVt = c{jj - Vt)dt + adWt, (51) 

which has the solution 

Vt = VqC-"' + ^(1 - e-"') + ae-"' f e^'dW,. (52) 

Jo 

The derivatives needed by the numerical methods are given in Table IXVI 



Model 


dVt 

at 


avt 

dWt 


Vasicek 


c(/x - Vt) 


a 


CEV 


k{9 - Vt) - \a^Vt 


aVt 



TABLE XV: equation array for Vasicek and CEV Models 

We set the parameters to c = .05, fi = .09, a = .03, and Vq = .08. We set the time step 
to dt = 2.4 and integrated to 12000. This is of course a very long dynamics. We plot the 
average relative error in Fig. 15 (a) for ANISE (solid curve) and SDE9 (dot-dashed curve). 
Both ANISE and SDE9 return results consistent with the requested tolerance. SDE9 returns 
a greater relative tolerance than that requested. 

The cpu times are compared in Table IXVll Here we see that SDE9 also runs somewhat 
faster than ANISE for this problem. 
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Model 


ANISE 


SDE9 


Ratio SDE9/ ANISE 


Vasicek 
CEV 


.79E+05 
.37E+05 


.57E+05 
.65E+05 


0.72 
1.75 



TABLE XVI: CPU times for tol = lO'^^ and 1 million trajectories. 
B. Mean- reverting CEV model 

Here the SODE is of the formal 

dVt = k{9 - Vt)dt + aVtdWt (53) 

which has the exact solution 

1 /■* 1 

Vt = exp{-{K + -a^)t + aWt] \Vq + kO ds exp{(/t + -a^)s - aWs}]. (54) 

2 JO 2 

The derivatives needed by the numerical methods are given in Table IX VI 

The parameters were chosen as k = .05, 9 = .09, a = .1, and Vq = .08. We set the time 
step to dt = 5 and integrated to 10000. The average relative error is plotted in Fig. 15 
(b) for ANISE (solid curve) and SDE9 (dot-dashed curve). Both ANISE and SDE9 return 
results consistent with the requested tolerance. Once again SDE9 returns a better relative 
tolerance than that requested. 

Table IXVII contains the cpu times for the two methods. SDE9 takes 75 % longer than 
ANISE. 

IV. CONCLUSIONS 

Good convergence is obtainable using both ANISE and SDE9 for all the problems con- 
sidered. In most cases ANISE runs roughly twice as fast. For the Vasicek model in section 
III SDE9 performed 40 % faster than ANISE. ANISE performed several hundreds of times 
faster than SDE9 for the Affine model in section II. 

In addition to our study of convergence, we examined the accuracy of individual tra- 
jectories for a given requested relative accuracy. We found that both methods returned 
trajectories with relative accuracies consistent with the accuracy requested, even for very 
long integration times. 
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F2_ave.a' i 'F3_ave.s 

'F2 ave.s' 'F3 ave.a' 

1 1 1 1 1 1 .-(4 4 I 1 1 1 1 1 1 1 1 1 1 

2000 4000 6000 8000 1 0000 1 2000 1 000 2000 3000 4000 5000 6000 7000 8000 9000 1 0000 



(a) ANISE (solid) and SDE9 (dot-dashed) (b) ANISE (solid) and SDE9 (dot-dashed) 

for Vasicek. for CEV. 

FIG. 15: Average error in volatility Vt for individual trajectories. 



Both algorithms appear to be sufficiently accurate for the models considered. ANISE 
performed better overall. The two methods appear capable of handling larger systems of 
equations with more Wiener processes, and could therefore prove to be valuable computa- 
tional tools for further research in finance. 
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